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In the context of Gross-Pitaevskii theory, we investigate the unconventional Bose-Einstein con¬ 
densations in the two-species mixture with p-wave symmetry in the second band of a bipartite 
optical lattice. A new imaginary-time propagation method is developed to numerically determine 
the p-orbital condensation. Different from the single-species case, the two-species boson mixture 
exhibits two non-equivalent complex condensates in the intraspecies-interaction-dominating regime, 
exhibiting the vortex-anti vortex lattice conhguration in the charge and spin channels, respectively. 
When the interspecies interaction is tuned across the SU(2) invariant point, the system undergoes 
a quantum phase transition toward a checkerboard-like spin density wave state with a real-valued 
condensate wavefunction. The influence of lattice asymmetry on the quantum phase transition is 
addressed. Finally, we present a phase-sensitive measurement scheme for experimentally detecting 
the UBEC in our model. 

PACS numbers: 03.75.Nt, 03.75.Lm, 05.30.Jp, 05.30.Rt 


I. INTRODUCTION 

Unconventional condensate wavefunctions of paired 
fermions are identified by nontrivial representations of ro¬ 
tational symmetry, in contrast to the conventional coun¬ 
terpart with vanishing relative orbital angular momen¬ 
tum (0AM). Exploration of unconventional condensates 
dates back to the investigations of the A- and 5-phases of 
the superfluid ^He [1-4] and later the spin-triplet pairing 
in Sr 2 Ru 04 [5-8], which are characterized by the forma¬ 
tion of Cooper pairs with 0AM of L = 1 and spin-triplet 
of 5 = 1. High Tc cuprates are another celebrated exam¬ 
ple whose pairing symmetry is [9, 10]. 

Recently, considerable discoveries, both theoreti¬ 
cal [11-25] and experimental [26-33], were reported on 
the single-boson condensation in the metastable high or¬ 
bital bands of an optical lattice. The wavefunctions of 
this archetype of unconventional Bose-Einstein condensa¬ 
tion (UBEC) are identified by the nontrivial representa¬ 
tions of the lattice symmetry group, which oversteps the 
physical scenario set by “no-node” theorem - an underly¬ 
ing principle of low-temperature physics stating that the 
many-body ground-state wavefunctions of Bose systems, 
including the superfluid, Mott-insulating and supersolid 
states, are necessarily positive-definite under general cir¬ 
cumstances [18, 34]. In consequence, the wavefunctions 
of UBECs can be rendered complex-valued, and thus 
spontaneously break the time-reversal (TR) symmetry 
[18], which constitutes a remarkable feature of UBEC. It 
is anticipated that UBECs can sustain exotic phenom¬ 
ena not seen in conventional BECs, such as the non¬ 
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trivial ordering of 0AM moment, BECs with nonzero 
momentum, half-quantum vortex and the spin texture 
of skyrmions [18]. It is also worth mentioning that the 
0AM moment formation still survives when system en¬ 
ters the Mott-insulating regime wherein the global U(l) 
phase coherence of superfluidity is no longer retained [19]. 

The experimental realization of single-species BECs in 
the second band, where the condensed atoms survive a 
long lifetime before tunnel to the nearly empty lowest 
band [28-32] , has marked an important progress towards 
the creation and manipulation of UBECs in ultracold 
atoms. Depending on the lattice asymmetry, the time- 
of-flight (TOE) measurement revealed signatures of both 
real and complex condensates with p-wave symmetry and 
a large scale spatial coherence. The complex wavefunc¬ 
tions exhibit the configuration of a vortex-anti vortex lat¬ 
tice with nodal points at vortex cores as theoretically pre¬ 
dicted. More recently, a matter-wave interference tech¬ 
nique was employed to provide direct observations of the 
phase information of the condensate and to identify the 
spatial geometry of certain low energy excitations [32]. 
The realization of UBECs in even higher bands was also 
reported [29, 30]. 

In this work, we present a theoretical study of the 
UBEC in a two-species boson mixture where both species 
are equally populated in the second band of a bipartite 
optical lattice [28] . Our study initiates the search of new 
types of UBECs enriched by coupling spin degrees of free¬ 
dom with U(l) symmetry, TR symmetry and nontrivial 
representations of the lattice symmetry groups. To de¬ 
termine the wavefunction of the UBEC in the context of 
Gross-Pitaevskii equation (GPE), we develop a numeri¬ 
cal scheme which resorts to precluding the s-orbital com¬ 
ponents from the condensate wavefunction during the 
imaginary-time evolution of the full Hamiltonian. This 
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scheme enables us to determine the phase diagram of 
UBEC in a wide range of parameters corresponding to 
the inter- and intraspecies interaction. We find that the 
emergent phases of UBEC involve the Px ^ipy (complex¬ 
valued) and Px ± Py (real-valued) types of orbital order, 
which appear in different regimes of interaction that can 
be described as a consequence of the interplay between 
0AM and interaction energies, as will be discussed later. 

This paper is organized as follows. The Sec. II, we 
briefly account for the experimental setup of the bipar¬ 
tite two-dimensional lattice potential used in our model, 
including the symmetry analyses of the lattice configu¬ 
ration. The structure of the single-particle dispersion of 
the p-band is demonstrated. In Sec. Ill, the numerical 
implementation of the modified imaginary-time propaga¬ 
tion method is described, which, together with the Bloch 
wave approximation, enables to solve the GPEs in high 
bands. In Sec. IV, we explore the properties of UBECs 
and phase transitions in the symmetric and asymmetric 
lattices. Einally, a scheme for experimentally exploring 
the formation of UBECs in our model is addressed in 
Sec. V and conclusions are made in Sec. VI. 


II. THE OPTICAL LATTICE AND BAND 
SPECTRUM 


We consider the two-species BEC in the first excited 
orbital band of the bipartite optical lattice employed in 
the experiments [28], where the unit cell consists of two 
sites with different potential depths. The optical poten¬ 
tial V{x,y) is described by 


V{t) = —^ \r] [{bz cos a -h Bz sin a) 


( 1 ) 


where the unit vectors and By constitutes the basis 
of the light polarization; Vq is determined by the laser 
power; ki = 27 t ja^ is the laser wavevector; a is the po¬ 
larization angle with respect to 2 ;-direction; e is the re¬ 
flection loss; the intensity and phase differences between 
laser beams along the x- and ^-directions are described 
by T] and respectively. The symmetry analysis of the 
lattice configuration and the subsequent band structure 
calculations have already been presented in Ref. [19]. Be¬ 
low we recap this analysis in detail to make the paper 
self-contained. 

Eor the ideal case with 77 = 1, e = 1, and a = 0, the 
lattice potential is simplified as 


U(r) = —Vq (^cos‘^ kix ^ cos‘^ kiy 

+2 cos kix cos kiy cos 0 


tice. When 77 < 1 and e = 1, the lattice potential becomes 
V (r) = — Vb cos^ kix + cos^ kiy 

-\-2r] cos kix cos kiy cos 0^ , 

which still possesses the reflection symmetries with re¬ 
spect to both X and 7/-axes, but the point group symme¬ 
try is reduced to the orthorhombic one. Eor the realistic 
case with 77 < 1 and e < 1, the orthorhombic symmetry is 
broken such that in general no special point group sym¬ 
metry survives. Nevertheless, the lattice asymmetry can 
be partially restored at = cos“^ e, where the lattice 
potential becomes 


U(r) 


-^{(1 + r 7 ^)(l + 6 ^) + 2e^rf cos2kix 
+ 2 e^ 77 ^ cos2kiy 


+ 4 e 77 cos 2kix 


e cos{kiy - 6 >) + cos{kiy + 6 >) | 


and the reflection symmetry with respect to the 7/-axis is 
retrieved. Therefore we call the case of a with a ^ ao 
as “asymmetric” and that oi a = as “symmetric”, 
respectively. The lattice structure with the experimental 
parameters Vb = 6.2E^, 77 = 0.95, 0 = 95.4°, e = 0.81, 
and Of = 7 r /5 is shown in Eig. 1 (a). 
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FIG. 1. (a) The double-well optical lattice with the experi¬ 
ment parameters Vb = 6.2Er, 77 = 0.95, 0 = 95.4°, e = 0.81, 
and a = 7 r/ 5 . The white dashed line illustrates the half wave¬ 
length of laser ao and is the lattice constant. The A 

and B sublattice sites are denoted in (a), (b) The energy 

spectra of the second band, whose energy minima are located 
at Ki = ( 7 r/ 2 ao, 7 r/ 2 ao) and K 2 = (— 7 r/ 2 ao, 7 r/ 2 ao). (c) and 
(d) are the density profiles (upper panel) and phase profiles 
(bottom panel) for non-interacting gas for Ki and K 2 , re¬ 
spectively. 


which possesses the tetragonal symmetry. Since 0 con¬ 
trols the relative depth of the double-well inside the unit 
cell, tuning 0 away from 90° results in the bipartite lat- 


The Bloch-wave band structure of the Hamiltonian 
Ho = —h^V^/ 2 M + U(i*) can be calculated based on 
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the plane-waves basis. The reciprocal lattice vectors are 
defined as Gm^n = mbi +nb 2 with bi ^2 = (±7r/ao, 7r/ao) 
where ao is a half wavelength of the laser. The diagonal 
matrix elements are 


(k + Gm,n\Ho\]<. + Gm,n) — 

dokx 


X 


+ (m — n) 


+ 


CLoky 


+ (m + n) 


( 2 ) 


where k is the quasi-momentum in the first Brillouin 
zone, and the off-diagonal matrix elements are 

{kV|k + G:pi,o) = + cos 

(k|y|k + Go,±i) = + cosae"'"*'’), 

( 3 ) 

{kV|k + G;pi,±i) = -^r/'^ecosa, 

(k|Vk + G^i,^i)=-^e. 

The energy spectrum of the second band of the optical 
lattice, Eq. ( 1 ), is shown in Fig. 1 {b). Several obser¬ 
vations are in order. Firstly, The energy minima are 
located at Ki ^2 = bi, 2/2 with the corresponding wave- 
functions t/^Ki and ' 0 K 2 - For the symmetric lattice, t/^Ki 
and are degenerate due to reflection symmetry, while 
for the asymmetric lattice, the degeneracy is lifted. Sec¬ 
ondly, there are four points in the Brillouin zone (BZ), 
namely, the zero center O, the high symmetry point X, 
(7r/ao,7r/ao), and Ki^ 2 , which are TR invariant because 
their opposite wavevectors are equivalent to themselves 
up to reciprocal lattice vectors. As a result, their Bloch 
wavefunctions are real, in other words, they are stand¬ 
ing waves instead of propagating waves. Thirdly, the 
hybridized nature of t/^Ki and is also manifest in 
real-space: their wavefunctions are mostly in the super¬ 
position of the local s-orbital of the shallow well and the 
p-orbital of the deep well, which possess nodal lines pass¬ 
ing through the centers of the deeper wells as shown in 
Fig. 1 (c) and (d). 


III. THE MODIFIED IMAGINARY-TIME 
PROPAGATION METHOD 

In current experiments [28], the correlation effects 
are relatively weak due to the shallow optical potential 
depth, and thus the two-species UBEC can be well de¬ 
scribed by the coupled GPE as 


E^p{r) 


H°p+ 9pa\^c.{r)f 

a=A,B 


^/3(r) 


( 4 ) 


where = {-h‘^V^)/2Mp + l/(r) is the one-particle 
Hamiltonian and the wavefunction T /3 is normalized to 


the area of one unit cell, |T^(r)p = Q = 2 ao; §a /3 = 

9 a^ni 3 with n /3 the particle number per unit cell and gai 3 
the interaction strength between a and (3 species. 

Furthermore, in terms of and the real- 

space spin density distribution is defined as S(r) = 
(l/ 2 )T’t'(r)(jT(r) where T = and a de¬ 

notes the Pauli matrices in vector form. Explicitly, 
the Cartesian components of the spin density are re¬ 
lated to and by 3x3- iSy = and 

Sz = Obviously, the orientation of 

spin in xy plane depends only on the global phases of 
and 

In solving Eq. (4), we assume qaa = gsB, 9 ab = Qba, 
Ma = Mb = M [35, 36], and ua = tlb- Since the band 
minima are located at Ki ^2 5 we expand the two-species 
condensate wavefunction in terms of t/^Ki and 

/ ^'^(r) \ ^ ( C0S(5A^/’Ki(r) + sin(5A?/)K2(r) \ 

V 4'B(r) J \ COS(5BV’Ki(r) + silK^B'i/jKa W J 

( 5 ) 

In general, '0 Ki(i*) and are determined by the 

renormalized lattice potential, and are thus different from 
those based on the free band Hamiltonian Hq [37]. Be¬ 
cause the particle number of each species is conserved 
separately, the formation of two-species BEC sponta¬ 
neously breaks the U(l)xU(l) symmetry, leaving the 
freedom of choosing the condensate wavefunction by indi¬ 
vidually fixing the phase factor of t/^Ki (i*) in each species 
of Eq. (5). 

The theoretical model in the single-species UBEC 
based on the GP description has been investigated with 
a self-consistent approach [19, 22]. For the two-species 
case, the structure of competing orders is even richer 
than that of the single-species case. In the enlarged phase 
space, the orbital states can entwine with spin degrees of 
freedom. We introduce a modified imaginary-time prop¬ 
agation method to solve the two-species UBEC, which 
liberates us from the restriction of certain types of solu¬ 
tions and can be generalized to other higher orbital bands 
as well. Since the ordinary imaginary-time propagation 
method only applies to yield the ground-state conden¬ 
sate, in order to reach the UBEC in the second band, the 
new method is devised to constantly project the lower or¬ 
bital components out of the evolving (in imaginary time) 
condensate wavefunction, forcing the initial wavefunction 
evolve to the stationary solution in the target orbital. We 
have examined this method for one- and two-dimensional 
harmonic oscillators, and the resultant wavefunctions not 
only converge to the exact solutions, but also yield the 
correct degeneracy of high energy levels. 

The implementation of our imaginary-time propaga¬ 
tion algorithm is summarized as follows. We start by 
initializing a trial condensate wavefunction in the form 
of Eq. (5) with ' 0 Ki ,2 determined by U(r) of the empty 
lattice. After the propagation of one time step, we ar¬ 
rive at a new set of and which are then em¬ 
ployed to generate the renormalized lattice potential 
Veff,a{^) = 4^(r) + E/?5/3a|^a(r)P- Then we solve the 
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s-orbital states |(/:’k)a at k = Ki and K2 based on Veff^ai 
and construct the projection operator 

P=^— |^k)aa(^k|- (6) 

a=A,B k=Ki,K2 

After projecting out the 5 -orbital component by apply¬ 
ing P to we proceed to the next step of imaginary- 
time evolution. The above process is repeated until the 
convergence is achieved and P is updated in each step. 
To assure the reliability of this method, we choose sev¬ 
eral different initial trial wavefunctions and add small 
complex random noises to break any specific symmetry 
which could lock the solution. Every simulation was im¬ 
plemented with a sufficiently long time to ensure that 
the energy converges. We have successfully reproduced 
the one-species UBEC solutions in the second band, and 
confirmed the results consistent with the previous works 
[19, 22]. The interaction strengths are much smaller than 
the energy difference between the s and p-orbital bands 
in our simulations, and thus the band mixing effect is 
negligible. 


IV. MAIN RESULTS 
A. The symmetric lattice 

We first consider the symmetric lattice and the compe¬ 
tition between intra- and interspecies interactions which 
determines the condensate wavefunctions. Defining 7 = 
QabIqaa-, we start with an SU( 2 ) symmetry breaking 
case in the regime of 7 < 1. When 7 = 0 , the system 
simply reduces to two decoupled single-species problems 
and each of them is in the complex condensate exhibit¬ 
ing nodal points rather than nodal lines. Accordingly, 
there are two p-orbital condensations characterized by 
substituting the following phase angles into Eq.(5): (I) 
, Sa = ^B = f, 

/ a J_ / V’Ki(r) +#K2(r) \ 

V^B(r)y ^/2 V V’Ki(r)+#K2(r) y ’ 

and (II) (pA = —= ±f, Sa = Sb = f, 

^A(r) \ J_ / V'Ki(r) +#K2(r) \ 

J V2 V V’Ki(r) -#K2(r) / ■ 

When 0 < 7 < 1, the corresponding p-orbital solutions 
take the forms of state (I) and (II) as well. 

States (I) and (II) possess different symmetry struc¬ 
tures, as illustrated in Eig. 2 (a) to (d). Species A and 
B can be interpreted as a Kramers doublet, and a com¬ 
monly used Kramers-type TR transformation is defined 
as T = io-yC where C is complex conjugation opera¬ 
tion and ay is the Pauli matrix. T keeps particle num¬ 
ber and spin-current invariant but flips the sign of spin 
and charge current, and it satisfies = —1. Eor state 


(I), its axial 0AM moments of two species are paral¬ 
lel exhibiting a vortex-antivortex lattice configuration, 
the condensate spin is polarized along the x-direction, 
which obviously breaks Kramers TR symmetry. As for 
state (II), its axial 0AM moments are antiparallel to each 
other exhibiting a spin-current vortex-anti vortex lattice 
configuration. Although spin current is invariant under 
Kramers TR transformation, the spin density exhibits 
the in-plane spin texture with the winding number ±2 
around each vortex core, which also breaks the Kramers 
TR symmetry. Nevertheless, state (II) is invariant by the 
anti-linear transformation T' = which is equivalent 
to a combination of the TR transformation followed by a 
rotation around the z-ax.is at tt. Since T'’^ = 1 , it is no 
longer a Kramer transformation, which maintains the xy- 
components of spin invariant but flips the ^-component 
of spin. 

States (I) and (II) give rise to the same particle den¬ 
sity and kinetic energy distributions for both species, and 
thus their energy are degenerate at the mean-field GPE 
level. Nevertheless, since they are not directly connected 
by symmetry, this degeneracy is accidental and only valid 
at the GPE level. The system symmetry allows a current- 
current interaction between two species, which is absent 
in the bare Hamiltonian, but could be effectively gener¬ 
ated through quantum fluctuations for low energy physics 
in the sense of renormalization group. Since the current 
density distributions of two species are the same in one 
solution but are opposite in the other. This emergent in¬ 
teraction would lift this accidental degeneracy. However, 
this is a high order effect beyond the GPE level, which is 
certainly an interesting subject for future investigations. 

The spatial distributions of the population and phase 
of both condensate species, together with the correspond¬ 
ing spin texture are shown in Eig. 2 (a) and ( 6 ), respec¬ 
tively. The particle density mainly distributes in the 
shallow sites which is the nodeless region corresponding 
to the 5-orbital, while the density in the deep sites at 
which the nodal points are located corresponding to the 
Pa,-orbitals. Each species exhibits a vortex-antivortex 
lattice structure: The vortex cores are located at the deep 
sites, and the nodeless region exhibits the quadripar¬ 
tite sublattice structure featuring the cyclic phase factors 
exp(i 7 rn/ 2 ) for n G {1,2,3,4} in the shallow sites. Eor 
state (I), both species exhibit the same vorticity distribu¬ 
tion and thus the spin density orientation lies along the 
x-direction according to the phase convention of Eq. (5). 
There is no preferential direction of spin orientation in 
xy plane due to U(l) symmetry generated by the total 
2 )-component spin. Eor state (H), the two species ex¬ 
hibit opposite vorticities, and the configuration is a spin- 
vortex-antivortex lattice. In both cases, the vorticity or 
the spin vorticity patterns exhibit a double period of the 
lattice potential. 

With 7 = 1 , the sum of interaction energies is ren¬ 
dered an SU( 2 )-invariant form such that the wavefunc¬ 
tions of UBEC become highly degenerate. At this point, 
the states (I) and (H) persist as expected. Because of 
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FIG. 2. Spatial distribusions of the density, phase, and spin 
texture of the condensate wavefunctions corresponding to var¬ 
ious states in the symmetric lattice are showcased in different 
groups of subplots: state (I) [(a) and (b)], state (II) [(c) and 
(d)] and checkerboard stdite [(e) and (f)]. The intra-species in¬ 
teractions are qaa = Qbb = 0.025Er with Er = h^kf/2M; 
the interspecies ones are qab = 0.25^aa for (a) and (5), 
and QAB = for (c) and (d), respectively. In (a), (c) 

and (e), the density and phase distributions are shown in the 
upper and lower panels for each species, respectively. The 
spin texture configurations are shown in (6), (d) and (/), 
respectively, with arrows indicating the orientation of spins 
and color bars representing the values of Sz- The parame¬ 
ters used are Vo = 6.2Er, rj = 0.95, 0 — 95.4°, e = 0.81, 
a = (Ao = cos“^ e 35.9°. 


the SU( 2 ) invariance, we can further apply the global 
SU( 2 ) rotations to states (I) and (II) [38]. For state (I), 
the constraint of maintaining ua = tib does not allow 
new states under the form of Eq. (5). For state (II), 
any SU( 2 ) rotation still maintains ua = tlb- For ex¬ 
ample, after a rotation of — 7 r /2 around the ^-axis, we 
arrive at = (V^Ki, —arid a subsequent 

7r/2-rotation around the x-axis yields 

( ^Air)\ ^ J_ [ V’Ki(r) + \ /O'! 

) y2 V V’Ki(r) - / ■ 

Next we consider the case of 7 > 1, where the degen¬ 
eracy of the SU( 2 ) invariant condensate wavefunctions is 
lifted. In this case, within the convention of Eq. (5), the 


solution of Eq. (9) is selected, whose density, phase and 
spin distributions are plotted in Fig. 2 (c) and (d). We 
see that bosons of different species occupy mostly the 
shallow sites in a checkerboard pattern with staggered 
spin density distribution. The condensate wavefunction 
in each species becomes real-valued with square-shaped 
nodal lines along with the period-doubled density pro¬ 
file, and we call this configuration the checkerboard state. 
In the single-species case [19, 22], the real non-Bloch 
states t/^Ki (r)±' 0 K 2 (i*) are always more energetic than the 
complex non-Bloch states t/^Ki (r) ± '<^'0K2 (i*) and the real 
Bloch states t/^Ki and 7 /^X 2 ? because the density distribu¬ 
tions of the real non-Bloch states are less uniform than 
those of the latter ones. However, the conclusion is oppo¬ 
site in the two-species case: both species exhibit strong 
constructive and destructive interferences between t/^Ki 
and 7 /^X 2 alternatively in adjacent shallow sites, and their 
real-space density distributions avoid each other and ex¬ 
hibit the checkerboard pattern. Consequently, the dom¬ 
inant interspecies interaction is greatly suppressed and 
the checkerboard state turns out to be the least ener¬ 
getic. 

In the strongly repulsive regime (7 > 1), however, it 
is possible that the system develops isolated ’’ferromag¬ 
netic” single-species domains. The case of spatial sep¬ 
aration has been discussed for the bosonic mixture in 
the s-orbital bands of optical lattices in the same in¬ 
teraction regime [39-42]. When this scenario occurs in 
p-orbital bands, isolated domains of either species can 
choose themselves in whichever of the complex states, 
' 0 Xi(i*) ± i' 0 x 2 (i*)- We call such a configuration the 
spatially phase-separated spin-polarized state. Seem¬ 
ingly, this state could have an energy lower than that 
of the checkerboard state because of the vanishing in¬ 
terspecies interaction. In Fig. 3, we plot E/Nfoi of 
full spin-polarized state with the complex condensate 
'0Ki(r) +i' 0 x 2 (^^): or, its TR breaking counterpart. Sim¬ 
ple numerical test shows that the energy per particle 
of the checkerboard state is very close to that of the 
fully spin-polarized state. When the initial state is 
prepared with ua = tib^ the fully spin-polarized state 
becomes phase-separated spin polarization accompanied 
with the formation of inhomogeneous ferromagnetic do¬ 
mains, which cost the domain energy. In spite of that, 
the checkerboard state of Eq. (9) is still the prevailing 
UBEC state in this regime. Another issue is the time 
scale: Starting from the unpolarized initial state, form¬ 
ing ferromagnetic domains is a process of phase separa¬ 
tion with a large scale arrangement of real space boson 
configurations. It is much longer than the time scale of 
the formation of the checkerboard state which only needs 
local phase adjustment of boson configurations. 


B. The asymmetric lattice 

Next we consider the interplay between lattice asym¬ 
metry and interactions. The lattice asymmetry breaks 
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FIG. 3. E/Ntot v.s. qabIqaa in the symmetric lattice with 
a = qq. Red dot, blue triangle, and green square are for 
different values of qaa = 0 . 015 F ^ r , 0 . 025 F ^ r , and 0 . 05 F ^ r , 

respectively. The dashed horizontal line represents the energy 
for 7 > 1 without i ncluding domain walls, z.e., solving the GP 
equation assuming fully polarization. The parameter values 
are the same as those in Fig. 2. 



FIG. 4. The condensate fraction of '0Ki, cos^ S, as a function 
of 7 = QabIqaa- The red dot, blue triangle and green square 
indicate qaaIEt — 0 . 015 , 0.025 and 0.05 respectively, and the 
dotted lines depict the fraction of spin-polarized state with 
the same total particle numbers for each interaction strength. 
The parameter values of the optical lattice here are the same 
as for Fig. 2 except for a = 7r/5 > ao- 


the degeneracy between the single particle states t/^Ki 
and ' 0 K 2 - Without loss of generality, we choose a > 
which sets the energy of slightly lower than that of 
t/^Ki , such that the calculated condensate wavefunctions 
satisfy 5a = ^ ^In Fig. 4, the condensate 

fraction of t/^Ki , cos^ 5^ is plotted as a function of 7 at 
various values of Qaa- We find that the lattice asymme¬ 
try effect is more prominent for weak interactions. At 
9aa = tfie condensate fraction of ' 0 Ki(i*) van¬ 

ishes when 7 < 0.5. The corresponding density, phase, 
and spin density distributions are depicted in Fig. 5 (a) 
and ( 6 ). This is a real Bloch-type UBEC with a stripe¬ 
like configuration and an in-plane spin orientation. With 
increasing 7 , t/^Ki and 7 /^X 2 superpose in a complex way 
with (pA = 4^3 = ± 7 r /2 or 4>a = —4b = i 7 r/ 2 , but 



FIG. 5. (a) Density, phase and (b) spin texture of qaa = 
Qbb = 0 . 015 Er , and 7 = 0 . 1 . The parameter values of the 
optical lattice here are the same as for Fig. 1 in the main text 
except a = 7 r/ 5 . 


cos^ 5 remains small even at 7 = 1. We note that, only 
when 7 > 1 , does the condensate quickly evolve to the 
checkerboard state. As qaa increases, the complex non- 
Bloch condensates become more and more prominent, as 
shown in Fig. 4. 

Hitherto, we have concluded that the two-species p- 
orbital condensation can manifest itself in different forms 
of non-Bloch condensation: the complex states (I) and 
(II), the real checkerboard state in Eq. (9) and the spa¬ 
tially phase-separated spin-polarized state. Since all 
these states are linear combination of ipKi (r) and 7 /^X 2 (i*)^ 
it is expected that four Bragg maxima would develop 
around the quasi-momenta, ±Ki ^2 in the TOE spec¬ 
tra [12, 19, 28]. Given the condensate fractions 4'Ki of 
Eig. 4, the states (I) and (II) as well as the spatially 
phase-separated complex spin-polarized state show that 
the relative intensities of these two pairs of peaks are 
dependent on the lattice asymmetry. However, when 
7 > 1 , the condensate fractions of t/^Ki and 7/^X2 foi* fhe 
real checkerboard state quickly become nearly equally 
populated and thus the Bragg peaks of the TOE spec¬ 
tra have almost equal intensities, irrespective of the lat¬ 
tice asymmetry. This experimental observation could di¬ 
rectly exclude the phase-separated spin-polarized state 
and provide supporting evidence for the phase transition 
from the complex UBECs towards the real-valued UBEC 
driven by the interspecies interaction. 


V. EXPERIMENTAL SCHEME FOR PHASE 
MEASUREMENT 

The two-species UBEC can be realized and observed 
by state-of-art experimental techniques [28-32]. Utilizing 
two different hyperfine spin states of an atom (labeled as 
the A- and H-species) [35, 36], one first creates a conden¬ 
sate of sole species in the superposition of ' 0 Xi(r) and 
' 0 K 2 (i*) which are the degenerate lowest-energy states in 
the p-orbital band. A 7 r/ 2 -Raman pulse is applied to con¬ 
vert half of the already condensed atoms into the other 
species. After tuning the interspecies atomic interaction 
with Eeshbach resonance, the system is held for some 
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time to let it relax to the intended non-Bloch p-orbital 
states, ^a,b = ('^Ki + whose phase infor¬ 

mation can be inferred by matter-wave interferometry as 
explained below. 

After the preparation of the two-speicies condensate, 
the atoms are then released from optical lattices and sub¬ 
sequently experience a Stern-Gerlach splitting during the 
ballistic expansion. Precisely, by applying a pulsed mag¬ 
netic field gradient, the atoms are accelerated by a spin- 
dependent force [43], F /3 ocm/ 3 | 5 |z (m /3 is the projection 
of spin), and thus the two-species UBEC breaks into spa¬ 
tially separated parts along 2 : direction. A second 7 r /2 
Raman pulse is then applied to mix states of different 
momenta, leading to 



where 4> accounts for the accumulated phases for the dy¬ 
namical effects involved, and Pa,s denote the momenta 
acquired by atoms after the Stern-Gerlach splitting. At 
this stage, the motion of each species is described by 
a wavepacket consisting of a superposition of two non- 
Bloch states with different quasi-momenta which inter¬ 
fere with each other along z direction during the TOE 
[32, 43]. The phase difference = 4>a — 4^b can be 

inferred from the interference patterns imaged along ver¬ 
tical and horizontal directions for each species, as demon¬ 
strated in [32]. It can be shown from the Eq. (10), that 
among the Bragg maxima, the Ki and K 2 columns pos¬ 
sess the same interference pattern, except the positions 
of fringes in the two columns are shifted by a phase an¬ 
gle, \/S.(j)AB\‘ By comparing the positions of fringes in the 
Bragg peaks, one can expect, when 7 < 1 , = 0 

for state (I) and \/S.(j)AB\ = tt for state (II). Our scheme 
provides a feasible way for phase measurement in the 
current system. 

VI. SUMMARY AND DISCUSSIONS 

In summary, we have studied the two-species p-orbital 
BECs in the experimentally accessible regime by a new 
imaginary-time propagation method for coupled GPEs, 
which can be applied to solve UBECs in higher bands. 
The competition between inter- and intraspecies inter¬ 
actions drives the transition from two non-equivalent 


complex-valued states, possessing respectively broken 
and unbroken TR symmetry, to a real-valued checker¬ 
board state with a staggered spin density structure. We 
have also proposed experimental schemes to study the 
UBECs of the mixture. The current study paves the way 
for approaching the least explored p-orbital physics of 
multi-species bosonic systems. Our theory can be also 
generalized to study the superfluity and magnetism of 
spinful p-orbital condensation in the presence of spin- 
dependent optical lattices or exotic spin-exchange inter¬ 
actions. 

We have used the GPE method throughout this arti¬ 
cle, whose applicability is justified in the limit of weak 
inter-species interaction. In this case, the two-species 
problem studied here is reduced to two weakly cou¬ 
pled single-species problems, for which previous works 
show that the GPE equation has captured the essen¬ 
tial physics of the complex p-orbital condensates being 
the energy-minima. When the interspecies interaction 
becomes stronger, however, the entanglement between 
two-species would become important. In this case, in¬ 
deed more exotic states beyond the GPE level is also 
possible. Eor example, the singlet paired boson con¬ 
densation, whose spatial pair-wavefunctions are antisym¬ 
metrized, and thus reduces the inter-species repulsion. 
This state is highly entangled and beyond the GPE equa¬ 
tion level. Nevertheless, the mean-field GPE equation is 
still a natural beginning point on this challenging prob¬ 
lem. The checkerboard state already investigated in this 
article remains a potential competing state, then both 
species avoid each other in their real-space density distri¬ 
butions characterized by a staggered spin-density struc¬ 
ture, which also greatly reduces the interspecies repul¬ 
sion. We would leave a detailed study on novel states 
beyond the mean-field GPE level and their competitions 
with the single-boson condensate for a future publication. 
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